Synthetic

We create a “synthetic population” by matching individuals surveyed for travel to the individuals surveyed for physical activity, based on age and gender. Then each person in the synthetic population has a set of trips taken and an amount of non-travel physical activity. This population is assumed to be representative of the population of the city or setting under study.

Scenarios

From the synthetic population we create “scenarios” of different pictures of travel by changing certain trips (e.g. swapping the mode of travel for particular trips and/or particular people). For the purposes of illustration, we choose a simple scenario for Accra in which every person gains one one-kilometre walk.

Distances

For the ITHIM, we process the travel data into distance data of various formats. First, we augment PT trips with a walk component: we take the set \(\{T_j:m(j)=\text{bus}\}\), which has \(J\) entries. We add \(J\) journeys with mode “walk” and duration \(\epsilon\) to the set \(T\). (This step would not be necessary for a travel survey which already includes all stages for each trip.)

Then we summarise the total distances and durations. The distance set is labelled \(\hat{T}_{m,s}\), representing the total distance travelled per mode per scenario, and the duration set is constructed analogously. We make an assumption that the distance for bus drivers scales linearly with bus passengers, where the distance in the baseline is defined based on the bus driver distance relative to car.

Distance and duration summaries include all person travel in the synthetic population and, in addition, some trips that are attributed to no-one in the synthetic population. These trips are used to represent “fleet” travel, which is uncaptured in travel surveys, and unchanged in scenarios. It’s important to include these trips because for modes such as motorcycle, where a considerable component of travel comes from fleet, attributing all the harms associated with the mode to only those surveyed dramatically overestimates the contribution of the individuals. In scenarios, then, we would be predicting emissions and injury rates that are far too high per person.

Background PM2.5

We use the total distance by mode (\(\hat{T}_{m,s}\)) to calculate the total PM2.5 in the scenario (\(\bar{P}_s\)). First we calculate the emission inventory in terms of fractions. If the confidence \(P\) is 1 then we use \(\dot{P}_{m}=\gamma_m\) directly. If \(P<1\), we sample the emission inventory fractions as follows: \[\tilde{\gamma}_m \sim \text{Dir}(f(P,\gamma_m)),\] \[\dot{P}_m=\tilde{\gamma}_m\] where we choose some mapping function such as \(f(P,\gamma_m)=\frac{\gamma_m}{\sum_m\gamma_m}10^{5P}\) (Figure 5).

Then we multiply vehicle distance by vehicle emission factor (where the emission factor is the emission inventory (fraction) divided by distance at baseline). \[\tilde{P}_{m,s}=\dot{P}_m\frac{\hat{T}_{m,s}}{\hat{T}_{m,s=\text{baseline}}}.\]

For the baseline, \(\tilde{P}_{m,s}\) is equal to the emission inventory fraction. For each scenario, it is scaled. We take the sum over modes to get a scalar for emissions in scenarios. \[\hat{P}_{s}=\sum_{m} \tilde{P}_{m,s}.\] Then the background PM2.5 in each scenario is the sum of the transport component and the non-transport component: \[\bar{P}_s=\eta(\zeta\hat{P}_s + 1-\zeta).\]

PM2.5 per person

Dose–response relationships are based on background concentrations in cities, and not on personal exposures. In our person-level model, we expect that individuals will have different exposures based on their different travel behaviours. While we can neither estimate individual exposures nor map them via any dose–response relationship to a relative risk, we aim to capture interpersonal differences via assumptions about the mechanistic relationships between behaviours and exposures, and anchor those to the background air pollution in their locality, for which dose–response estimates of relative risk exist.

Therefore, when we write “individual exposure,” it should not be read as an aetiologically relevant exposure, but rather it should be read as a unit- and scale-agnostic metric that maps an individual to a particular dose with reference to city-level measurements of air pollution. New research is required to elicit dose–response relationships between individual health exposures and health outcomes \[citation from JPO?\]. Our approximations act as assumptions for now and placeholders until more detailed data are available. \[Perhaps we should use words other than "individual exposure". How about "individual PM2.5 reference values"?\]

Individual exposures to PM2.5 are calculated using the background PM2.5 (\(\bar{P}_s\)) and the trip sets. There are three major components to daily exposure: one, a person’s total inhalation off road; two, a person’s inhalation on road in a vehicle, and, three, a person’s inhalation on road while undertaking active transport. Each category has an amount of background PM2.5 and a ventilation rate which together inform overall exposure.

The ratio of exposure off road to that on road is a function of total PM2.5, defined as \[K_s=C_4-C_5\log(\bar{P_s})\] as in Goel et al. (2015). This defines the exposure of a person in an open vehicle (i.e. pedestrian, cyclist or motorist): \[\tilde{K}_{m\in\{\text{walk,cyc.,mot.}\},s}=K_s\] and it is used to calculate in-vehicle exposure, assuming an exposure of \(C_2\) with the window closed, and a proportion \(C_3\) of vehicles having closed windows: \[\tilde{K}_{m\not\in\{\text{walk,cyc.,mot.,sub.}\},s}=C_2C_3+K_s(1-C_3).\] The exposure in a subway is constant and not dependent on the road ratio \(K_s\): \[\tilde{K}_{m=\text{sub.},s}=C_6.\]

Ventilation rates are calculated for each mode, assuming a base-level inhalation rate of \(C_1\): \[V_m=C_1+\frac{1}{2}\lambda_m,\] where \(\lambda_m\) is the MMETs for mode \(m\).¹ Then the air inhaled during travel per person is \[\tilde{V}_{i,m,s} = \smashoperator{\sum_{\substack{j:i(j)=i,\\m(j)=m,\\s(j)=s}}}T_j\cdot V_{m(j)}/Z_{m(j)}\] and the rate of inhalation during travel is \[\bar{V}_{i,s} = \sum_{m}\tilde{V}_{i,m,s}\cdot \tilde{K}_{m,s}.\] The air inhaled when not travelling is \[\hat{V}_{i,s} = C_1\left(24-\smashoperator{\sum_{\substack{j:i(j)=i,\\s(j)=s}}}T_j/Z_{m(j)}\right).\] Together, the PM2.5 exposure is calculated as the total PM2.5 inhaled per hour as follows: \[\check{W}_{i,s} = \frac{\bar{P}_s(\bar{V}_{i,s}+\hat{V}_{i,s})}{24}.\]

AP–disease dose–response relative risk

We use each person’s exposure to PM2.5 (\(\check{W}_{i,s}\)) to calculate their relative risk of five diseases (IHD, lung cancer, COPD, stroke, LRI), using curves parametrised by four disease-specific variables (Burnett et al. 2014). Of the five diseases, two (IHD and stroke) have parameters specific to age groups starting at age 25.² The other three (lung cancer, LRI and COPD) have one set of parameters for all ages.

The curves are in the form of samples of the set of four parameters. We model the densities of these samples (using a quantile for parameter 3, kernel density estimation for parameter 2, and GAMs for parameters 1 and 4)³ in order to draw either the median or random samples via their quantiles (\(\xi_{h,w}\)) as follows:

\[\begin{aligned} \tilde{H}_{a,h,w=3}=&\text{CDF}_{H_{a,h,w=3}}^{-1}(\xi_{h,w=3})\\ \tilde{H}_{a,h,w=2}=&\text{CDF}_{H_{a,h,w=2}|\tilde{H}_{a,h,w=3}}^{-1}(\xi_{h,w=2})\\ \tilde{H}_{a,h,w=1}=&\text{CDF}_{H_{a,h,w=1}|\tilde{H}_{a,h,w=2},\tilde{H}_{a,h,w=3}}^{-1}(\xi_{h,w=1})\\ \tilde{H}_{a,h,w=4}=&\text{CDF}_{H_{a,h,w=4}|\tilde{H}_{a,h,w=3},\tilde{H}_{a,h,w=2},\tilde{H}_{a,h,w=1}}^{-1}(\xi_{h,w=4}).\end{aligned}\]

From these parameters, the relative risk of mortality is defined \[\check{H}_{h,i,s} = 1 + \tilde{H}_{a=a(i),h,w=1}\times \left\{1 - \exp\left(-\tilde{H}_{a=a(i),h,w=2} (\check{W}_{i,s} - \tilde{H}_{a=a(i),h,w=4}) ^ {\tilde{H}_{a=a(i),h,w=3} }\right)\right\}.\] For \(h\not\in\{\text{IHD, lung cancer, COPD, stroke, LRI}\}\), \(\check{H}_{h,i,s} \equiv 1\).

Individual-level MMETs

Using trip sets and the synthetic population, we calculate total MMETs per person (\(M_{i,s}\)) as the sum of walking MMETs per day and cycling MMETs per day, which are scaled up to a week via the scalar \(\chi\), and work/leisure MMETs per week as follows: \[M_{i,s}=\smashoperator{\sum_{\substack{j:i(j)=i,\\s(j)=s}}} \chi\lambda_{m(j)}\frac{T_j}{Z_{m(j)}}+\theta Y_i.\]

Here, \(\lambda_m\) is the MMETs for mode \(m\), \(T_j\) the distance of trip \(j\), \(Z_m\) the speed of mode \(m\), \(\theta\) the scalar for background PA, and \(\check{Y}_i\) the amount of work and leisure (background) PA of person \(i\).

We calculate \(\check{Y}_i\) according to the confidence in the PA survey, \(\tilde{Y}\). First, we calculate the raw probability that a person in demographic group completes no non-travel PA: \[\dot{Y}_{a,g}=\text{Prob}(Y_i=0)|_{a(i)=a,g(i)=g}.\] Then, we set the probability to use, \(\bar{Y}_{a,g}\), as \(\dot{Y}_{a,g}\) if our confidence \(\tilde{Y}\) is 1. If \(\tilde{Y}<1\), we map the confidence to parametrise a Beta distribution somehow, e.g. (Figure 4) \[\begin{aligned} \hat{Y} &= 500^{\tilde{Y}+0.2},\\ \beta & = \hat{Y}\dot{Y}_{a,g}\left(\frac{1}{\dot{Y}_{a,g}} -1\right),\\ \alpha & = \hat{Y} - \beta,\\ \bar{Y}_{a,g} & = \text{CDF}_{\text{Beta}(\alpha,\beta)}^{-1}(\tilde{\theta}).\end{aligned}\]

Finally, for each person in the population, we sample non-travel MMETs as zero with probability \(\bar{Y}_{a,g}\) and from the raw non-zero density of their demographic group with probability \(1-\bar{Y}_{a,g}\).

PA–disease dose–response relative risk

We use each person’s MMETs per week (\({M}_{i,s}\)) to calculate their relative risk of six diseases, using curves found from meta analysis. Each disease (except type 2 diabetes) has a threshold beyond which there is no further change in relative risk. (Total cancer, breast cancer, colon & rectum cancer, endometrial cancer, & coronary heart disease: 35; lung cancer: 10; stroke: 32; all cause: 16.08.) For all individual cancers, we estimate the burden in terms of . For all other diseases, we use mortality.

The interpolation of the dose–response curve can be described as follows: \[\tilde{G}_{h,x}(m) = {G}_{h,x,z=\lfloor{m}\rfloor}+({G}_{h,x,z=\lceil{m}\rceil}-{G}_{h,x,z=\lfloor{m}\rfloor})\frac{m-\lfloor{m}\rfloor}{\lceil{m}\rceil-\lfloor{m}\rfloor}\] where we interpolate the mean, the lower bound, and the upper bound. These define the normal from which we sample, truncated at 0: \[F_S(v_x) = \mathcal{N}\left.\left(v_{x=2},\frac{v_{x=3}-v_{x=1}}{1.96}\right)\right|_{F_S(v_x)>0}.\]

Then the relative risk for each person for each disease is calculated for the quantile \(\phi_h\) as \[W_{h,i,s}=\text{CDF}_{F_S\left(\tilde{G}_{h,x}(M_{i,s})\right)}^{-1}\left(\phi_{h}\right).\] See Figure 2 for an example. For \(h\not\in\) {diabetes, total cancer, breast cancer, colon & rectum cancer, endometrial cancer, IHD, lung cancer, stroke, all cause}, \(W_{h,i,s} \equiv 1\).

\[Example of the dose--response workflow for PA and total cancer. Top left: results of the meta-analysis: mean, lower bound and upper bound of the relationship between PA dose and relative risk of disease. Top right: examples of 500 samples from the normal distribution defined by the mean of the top-left panel and the standard deviation defined by (upper bound minus lower bound) divided by 1.96 and truncated at 0. Bottom left: relative risks of individuals with PA doses 5, 25 and 45 are found by mapping their doses onto the median curve. This is used in the "constant" ITHIM use case. Bottom right: relative risks of individuals with PA doses 5, 25 and 45 for three different samples from the distribution shown in the top-right panel. For each sample, each individual is mapped onto the response curve it defines.\]

Processing

For the injury module, we require distance per travel mode per demographic group per scenario. It requires matching of mode names in the travel survey to mode names in the injury records, for example “taxi,” “shared auto,” “shared taxi” and “car” combine to form “car.” Some work is required to separate drivers from passengers; we currently assume all travellers are drivers, with the exception of “bus.”⁴

Modelling

We model injuries via regression by predicting the number of fatalities of each demographic group on each mode, which we calculate as a sum over all the ways in which they might have been injured (i.e., all the modes with which they might have collided).

The predictive covariates include the distances travelled by the parties and their demographic details. These requirements lead to a natural separation of the (training) dataset into two groups: a set for which we have distance data for the other party, and a set for which we do not have distance data for the other party. The former equates to a “who hit whom” (“whw”) matrix (albeit in a higher dimension), and will account for changes in injuries resulting of a change in strike-mode travel. The latter corresponds to causes of injury that will not change across scenarios, including “no other vehicle” and modes of transport that we do not consider to change, which might include trucks and buses if they are not somehow explicitly included in the trip set. We label this group “noov”: no or other vehicle.

Note on recategorization: for some of the modes (like for car, bus, motorcycle and cycle), if both strike and victim modes are same, we explicitly make the strike mode as ‘noov’

We model the number of injuries as a Poisson-distributed variable with an offset depending on the distance(s) and the reporting rate. We write the model as follows: \[\begin{aligned} \tilde{I}_{a,g,m_{\text{cas}},m_{\text{str}},s,y} \sim & \text{Poisson}(\bar{I}_{a,g,m_{\text{cas}},m_{\text{str}},s}); \\ % 'count~ \log\left(\bar{I}_{a,g,m_{\text{cas}},m_{\text{str}},s}|m_{\text{str}}\in\{m_{\text{whw}}\}\right) = & % cas_mode*strike_mode-log(injury_reporting_rate)+log(weight) {\kappa_0} +\sum_{i=1}^n\kappa_iX_i-\log(\sigma)+\nonumber\\ % +log(cas_distance) & \log(\bar{T}_{a,g,m_{\text{cas}},s}) % +(CAS_EXPONENT-1)*log(cas_distance_sum) +(E_{m_{\text{cas}}}-1)\log(\hat{T}_{m_{\text{cas}},s}) %+log(strike_distance)+(STR_EXPONENT-1)*log(strike_distance_sum) +E_{m_{\text{str}}}\log(\hat{T}_{m_{\text{str}},s});\\ \log\left(\bar{I}_{a,g,m_{\text{cas}},m_{\text{str}},s}|m_{\text{str}}\in\{m_{\text{noov}}\}\right) =& {\kappa'_0} +\sum_{i=1}^n\kappa'_iX'_i -\log(\sigma) % +log(cas_distance) + \log(\bar{T}_{a,g,m_{\text{cas}},s}) % +(CAS_EXPONENT-1)*log(cas_distance_sum) +(E_{m_{\text{cas}}}-1)\log(\hat{T}_{m_{\text{cas}},s}).\end{aligned}\] We choose this form of equation to enable (a) linearity in injuries with respect to total travel combined across modes and (b) linearity in injuries as subdivided by travellers of each mode.⁵ Note that this methodology is under development. See https://github.com/robj411/safety_in_numbers_power_law.

In general, one makes the best model that one can with the data available. We use, for Accra, the covariates casualty mode*strike mode, casualty age, and casualty gender, to form the model matrix \(X\). We also have another covariate, year, but our travel data are for one year only, so we model instead many observations of the same quantities.

We fit the coefficients \(\kappa\) of the model using data for ten years, assuming the same travel each year, which corresponds to the scenario \(s=\text{baseline}\). To incorporate the injury reporting rate (\(\sigma\)), we set it as an offset, whose value is 1 in fitting the model, and as specified in making predictions.

Prediction

We predict the number of injuries in each scenario based on the training model built from the baseline scenario. For Accra, we predict for the year 2016, using scenario-specific travel data. The number of fatalities is taken as the sum of fatalities over the strike modes: \[\check{I}_{a,g,m,o=\text{death},s}=\sum_{m_{\text{str}}}\tilde{I}_{a,g,m_{\text{cas}},m_{\text{str}},s}.\] We extrapolate injury deaths to injury YLLs via the ratio in the GBD: \[R_{a,g}=U_{a,g,h=\text{road injuries},o=\text{YLL}}/U_{a,g,h=\text{road injuries},o=\text{death}},\] \[\check{I}_{a,g,m,o=\text{YLL},s}=\check{I}_{a,g,m,o=\text{death},s}\cdot R_{a,g}.\]

PA and AP relative risks combined

We combine the relative risks of disease as a function of AP and disease as a function of PA through multiplication: \[\tilde{W}_{h,i,s}=W_{h,i,s}\cdot\check{H}_{h,i,s}.\] Not all diseases have a dose–response relationship for both AP and PA. Just stroke, lung cancer, and IHD have both. For the other diseases, only one of \(W_{h,i,s}\) and \(\check{H}_{h,i,s}\) is different from one.

Total disease burden

We calculate the total health burden relative to the reference scenario (which, for our Accra example, is the baseline) using the injury and health outputs combined with GBD data, via population attributable fractions (PAF). The PAF is defined: \[\hat{W}_{a,g,h,s}=\smashoperator{\sum_{\substack{i:a(i)=a,\\g(i)=g}}}\tilde{W}_{h,i,s}.\]

We then calculate the PAF relative to the reference scenario (“ref”) as: \[\bar{W}_{a,g,h,s}= \frac{\hat{W}_{a,g,h,s=\text{ref}}-\hat{W}_{a,g,h,s}}{\hat{W}_{a,g,h,s=\text{ref}}}.\]

We estimate the background burden of disease using GBD data and scaling based on the ratio of populations between country and city. For country populations \(N_{a,g}\), city populations \(\bar{N}_{a,g}\) and country burden \(U_{a,g,h,o}\), we estimate city burden as \[\bar{U}_{a,g,h,o}=U_{a,g,h,o}\frac{\bar{N}_{a,g}}{N_{a,g}}.\]

If we are scaling the background burden of non-communicable diseases (“Neoplasms,” “Ischemic heart disease,” “Tracheal, bronchus, and lung cancer,” “Breast cancer,” “Colon and rectum cancer,” “Uterine cancer”), we do so here. (Otherwise set \(\rho\equiv 1\).) \[\tilde{U}_{a,g,h,o}=\rho_h\cdot \bar{U}_{a,g,h,o}.\]

We combine the burden with the PAF: \[\hat{U}_{a,g,h,o,s}=\bar{W}_{a,g,h,s}\cdot\tilde{U}_{a,g,h,o}.\]

Injury burden

And, finally, for the injuries, we sum over modes to compute the burden, \[\bar{U}_{a,g,o,s}=\sum_m\check{I}_{a,g,m,o,s},\] and subtract from the values for the reference scenario: \[\hat{U}_{a,g,h=\text{road injury},o,s}=\bar{U}_{a,g,o,s=\text{ref}}-\bar{U}_{a,g,o,s}.\]

Physical activity

Each disease’s PA dose–response relationship is defined by a truncated normal. For each dose, there is a mean value, an upper bound, and a lower bound. For each person’s dose, we get the response by mapping the uniform random variable onto the truncated normal defined by the mean and bounds for that dosage. We use an R package called {drpa}, whose description can be seen here: https://meta-analyses.github.io/drpa/drpa.html

The relationships comes out of our work Garcia et al. (2018), where we look at the association of physical activity with various diseases.

Air pollution

For the AP relationship, there are four parameters per disease. We sample the first from an empirical distribution using the probability integral transform. We sample the second via the same method, conditioned on the value of the first, constructing their joint density with e.g. kde2d. The third parameter is sampled conditioned on the first and second, constructing their joint density using a GAM. The final parameter is sampled conditioned on the first, second, and third, constructing their joint density using a GAM.

As before, there is perfect correlation between individuals, i.e. if person A’s dose is greater than person B’s, then person A’s response is strictly greater than person B’s response.

The empirical distributions come from Burnett et al. (2014). There are four parameters per disease: IHD, lung cancer, COPD, stroke, and LRI. In addition, for stroke and IHD, there is a set of four parameters for each age group from 25 to 95 in five-year increments. In addition to our assumption that there is perfect correlation between individuals for diseases, we assume perfect correlation between ages for diseases. I.e., our four quantiles per disease will be applied to all age groups.

Summarise the travel survey

Our first step is to summarise the demographic information and the travel survey. In Table 5 are the demographic data: the populations of the city, and the number of people in the survey. In Table 6 are the summary statistics from the travel survey which we will use to generate synthetic travel. From the raw travel data (actually the adjusted trip set, which is the raw travel plus motorcycle, bus and truck trips), we calculate the probability to travel, \(B_{a,g,m}\), as the number of people who travelled by that mode divided by the number of people in the demographic group. It is to these values that we assign uncertainty through the variable PROPENSITY_TO_TRAVEL.

Adjust for trip weight

We define a trip-mode weight per mode, \(L_m\), so that a single trip by mode \(m\) in the baseline scenario represents \(L_m\) trips in reality. If in a scenario the mode of the trip is changed (\(m'\rightarrow m\)), it retains the weight of its original mode \(m'\). We use the trip-mode weight to (a) sample trips for scenario generation and (b) scale the probabilities to travel. Note that we don’t use it as trip weights are ordinarily used: to scale trip distance/duration. The result is a scaled set of probability parameters: \[\hat{B}_{a,g,m} = L_m{B}_{a,g,m}.\] (\(L_m\) should actually be the average trip-mode weight for all trips in the scenario. For baseline, it will be the same. It can/should differ for scenarios. Need to work out notation.)

Smooth probabilities

We estimate smoothed probabilities for any mode that has zero or one probability, using the rest of the probabilities for reference, replacing \(\hat{B}_{a,g,m}\). The reference to other modes is slightly problematic if e.g. in the scenario all travel increases. It means a 0 in the baseline is not the same as a 0 in the scenario. This could be fixed with more rigid coding.

Resample probabilities

We resample the adjusted probability parameters of Table 6. We describe each parameter with a Beta distribution as follows: \[\begin{aligned} \label{pitilde} \tilde{B}_{a,g,m} & \sim \text{Beta}(\alpha_{a,g,m},\beta_{a,g,m}); \\ \alpha_{a,g,m}+\beta_{a,g,m} & = 200; \\ \hat{B}_{a,g,m} & = \frac{1}{1+\frac{\beta_{a,g,m}}{\alpha_{a,g,m}}}.\end{aligned}\] Here, \(\tilde{B}_{a,g,m}\) is the new, variable parameter that we resample for each simulation. We sample via a quantile, \(\nu_m\sim\mathcal{U}(0,1)\), which we use to calculate \(\tilde{B}_{a,g,m}\). We set \(\alpha_{a,g,m}\) and \(\beta_{a,g,m}\) to sum to 200, making the distributions concentrated, and we constrain their relative values to enforce the distribution’s mean to equal \(\hat{B}_{a,g,m}\) (the right-most column in Table 6).

Inform individual travel

Each individual is assigned their own propensity to travel per mode, \(\delta_{i,m}\), which is a uniformly distributed random variable that is fixed and not resampled. Whether or not individual \(i\) travels by mode \(m\) depends on \(\tilde{B}_{a(i),g(i),m}\), which does vary: specifically, individual \(i\) travels by mode \(m\) if \(\delta_{i,m}<\tilde{B}_{a(i),g(i),m}\). (To mitigate the effects of these fixed random variables, a large population is required, e.g. more than 5000 individuals.)

\[Distributions for propensities to travel. The propensity is the probability that a member of a demographic group will travel by a particular mode on a particular day. In orange are the raw probabilities based on the synthetic trip set, here shown for the baseline scenario. Beta distributions are defined via a pointiness parameter, set to 200, such that $\alpha+\beta=200$.\]

(){#propensity_to_travel width=“70%”}